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Linear optical response 
of current-carrying molecular junction: 
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We propose a scheme for calculation of linear optical response of current-carrying molecular junc- 
tions for the case when electronic tunneling through the junction is much faster than characteristic 
time of external laser field. We discuss relationships between nonequilibrium Green function (NEGF) 
and time-dependent density functional theory (TDDFT) approaches, and derive expressions for op- 
tical response and linear polarizability within NEGF-TDDFT scheme. Corresponding results for 
isolated molecule, derived within TDDFT approach previously, are reproduced when coupling to 
contacts is neglected. 



I. INTRODUCTION 

Rapid progress of experimental capabilities in the field 
of molecular electronics necessiates development of the- 
oretical (and calculational) tools capable of explaining 
existing data and predicting (proposing) future experi- 
ments. While initially focus of both experimental and 
theoretical studies was on ballistic current-voltage char- 
acteristic of molecular junctions^ today more compli- 
cated phenomena are on the forefront of research. This 
includes inelastic transport (inelastic electron tunneling 
spectroscopy both far ofS and at resonance^), current- 
induced motionp nonlinear conductance (negative dif- 
ferential resistance* 5 - hysteresis* 6 - and switching 7 -), shot 
noise^ Coulomb blockade^ Kondo effect ) 10 ' 11 heating of 
molecular junctions ; 12 ! 13 etc. Recently optical properties 
of molecular junctions started to attract attention of re- 
searchers. Application of external laser field promises 
ability to effectively control transport properties of 
molecular devices^ while Raman spectroscopy, together 
with scanning tunneling microscopy (STM) and inelas- 
tic electron tunneling spectroscopy (IETS) can serve as a 
diagnostic tooh 15 ' 16 First experimental data in this direc- 
tion include light induced switching behavior of molecu- 
lar junctions ; 17 ' 18 ! 19 voltage effects on fluorescence^ - of 
molecules in nano junctions, and surface enhanced Ra- 
man scattering (SERS) from molecules positioned in nar- 
row gaps between metal nanoparticlesi 21 ' 22 ' 23 ' 24 Theo- 
retically opto-clcctronic properties of current carrying 
molecular junctions have been studied mostly within sim- 
ple models . 25 ' 26 ' 27 ' 28 Here we make a first step in direc- 
tion of ab initio calculations of optical properties of such 
junctions. 

Transport properties of molecular junctions are natu- 
rally described within non-equilibrium Green's function 
(NEGF) approach. 48 ' 49 A consistent way to treat trans- 
port of molecular junction within NEGF implies using 
many-body Green functions approach. However, com- 
plexity of the schemes limits applicability of this treat- 
ment to relatively simple molecular models only. 50 ' 51 ' 52 ' 53 
DFT is a well-established tool for accurate calculation 



of ground state electronic properties of isolated systems 
(atoms, molecules, solids) i 29 ' 30 ' 31 ' 32 DFT ability to deal 
with relatively extended realistic systems made the idea 
of merging the two approaches very appealing. First 
NEGF-DFT was implemented to steady-state ballistic 
(Landauer) transport calculations. 54 ' 55 ' 56 Later the cal- 
culations were extended to the case of inelastic elec- 
tron transport in the weak electron-phonon coupling 
limit. 57 ' 58 ' 59 ' 60 ' 61 For detailed discussion on different ap- 
proaches to transport through molecular junctions see 
Ref.Hl While NEGF-DFT calculations modedling IETS 
reported good agreement with experimental data ; 58 ' 59 re- 
sults for ballistic current-voltage calculations are not al- 
ways in quantitative agreement. Sometimes experimental 
data are reproduced by NEGF-DFT approach) 63 ' 64 while 
in other cases calculations yield a current that is orders 
of magnitude larger than experimental values. 65 ' 66 

Besides obvious uncertainties inherent for molecu- 
lar junction simulations, such as influence of contact 
geometr y 67 ' 68 and local environment 6 ^ on transport, a 
methodological DFT problem (when applied to trans- 
port situation) was pointed as a possible source of 
discrepancy. 70 ' 71 ' 72 ' 73 Possible errors can be grouped into 
several categories: 1. use of inappropriate exchange- 
correlation functional (ground state and/or spatially lo- 
cal character of functionals used, self-interaction errors, 
absence of derivative discontinuity and xc contribution to 
the electric field response); 2. inherent time-dependent 
character of transport, which goes beyond validity of 
DFT; 3. DFT (as well as TDDFT discussed below) is 
a theory which work for finite systems (external distur- 
bance local in space) only; 4. unphysical separation of the 
system into disconnected parts (contacts and molecule) 
at infinite past within NEGF. For a detailed discus- 
sion on application of DFT in transport calculations see 

Refs. [zHzizi. 

Optical response of system is defined by its electronic 
excitations. In order to calculate excited states two 
main approaches can be implemented. One is based 
on a many-body theory, where solution of the Bcthe- 
Salpctcr equation (BSE) is needed. Its computation- 
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ally expensive character limits applicability of the ap- 
proach to a relatively narrow range of problems (see e.g. 
Rcfs. 33.34.35). A much more numerically efficient al- 
ternative is the time dependent density-functional the- 
ory (TDDFT). For a review comparing these two tech- 
niques see Ref. [3(1 TDDFT is an extention of density 
functional theory ; 37 ! 38 which properly treats correlated 
excited states and where electronic excitations are asso- 
ciated with the poles of exact charge density response. 
TDDFT is known for successful calculations of optical 
spectra in many finite molecular system s 36 ! 39 Its numer- 
ical simplicity allows treating systems that involve hun- 
dreds of atoms, and a lot of work within the approach has 
been done for isolated molecules 40,41,42,43,44,45 ,46,47 Tak . 

ing into account size of realistic systems (molecules) used 
in molecular devices, TDDFT is the only tool available 
today capable of dealing with such calculations. 

Transport calculation schemes based on TDDFT were 
proposed as an alternative to NEGF-DFT method. Main 
differences are time-dependent character of the TDDFT 
scheme (TDDFT instead of DFT), absence of system 
partitioning, and finite size of the system under study. 
In particular, one of the proposed schemes, time depen- 
dent Kohn-Sham master equation approach, utilizes ring 
geometry, and linearly increasing magnetic field in the 
center of the ring provides driving force i 74 ' 75 Another 
TDDFT-NEGF approach considers time evolution of fi- 
nite linear system initially at equilibrium under influ- 
ence of external field (battery discharge) / 6 ! 77 ' 78 ! 79 While 
these approaches seem to solve (at least partially) prob- 
lems of NEGF-DFT, some questions still exist. For ex- 
ample, in the master equation approach, finite size of 
the ring forces to apply artificially large coupling be- 
tween electrons and the bath in order to achieve ther- 
malization in the part of the ring representing contact, 
which rises question of physicality of charge distribution 
in the contacts and, hence, current through the device 
part of the ring. TDDFT-NEGF also implements sys- 
tem of finite size, i.e. continuum character of states in 
the electrodes is questionable. Thus equivalence of tran- 
sient current calculated within the approach and realistic 
steady-state current is not obvious. The approaches are 
still have to be tested on the problems where NEGF- 
DFT failed. Note also that TDDFT has its own lim- 
itations (nonuniquencss of the excited-state potentials, 
question of stability and chaos of the mapping of densi- 
ties on potentials) i 80 ! 81 

While agreement on existence of methodological pit- 
falls of NEGF-DFT exists, the importance of those er- 
rors (i.e. if they are the cause of discrepancy between 
NEGF-DFT results and experimental data) is not com- 
pletely understood yetM Besides, within TDDFT-NEGF 
approach it was shown that after initial correlations die 
out and steady-state current is established, the last is 
given by Landauer-like formula with chemical potentials 
in Fermi distribution functions being shifted in accor- 
dance with extra exchange-correlation term originating 
from response to external field (bias) i 74 ' 76 The prob- 



lems of NEGF-DFT and TDDFT-NEGF schemes seem 
to stem from essentially different basic assumptions of the 
two (NEGF and DFT) theories, compatibility of those is 
not clear. However this is the only practical tool avail- 
able today capable of dealing with realistic simulations 
of molecular junctions. 

In this paper we consider optical response of current 
carrying molecular junction within NEGF-TDDFT ap- 
proach. We assume that initially steady-state current 
across the junction is established, and then the system 
is probed by an external laser field. The last is assumed 
to be a weak perturbation on top of the non-equilibrium 
steady-state. The assumption works in the case when 
timescale for electron transport is much shorter than 
the characteristic time of external field. In the opposite 
case (both times are comparable or the field is quicker) 
one would need to consider time-dependent transport, ei- 
ther within TDDFT-NEGFl 6 .! 7 !^!! 9 . or time-dependent 
NEGF— approach, explicitly. We postpone such consid- 
eration for future research. In treating steady-state flux 
we assume that Landauer formula is correct, while xc and 
chemical potential are adjusted properly, as is discussed 
in Refs. I74ll76l So, the treatment formally looks like the 
standard NEGF-DFT, however one has to keep in mind 
points mentioned above. 

In Section|TT]wc introduce model of molecular junction. 
Section IIIII briefly discusses general differences between 
NEGF and TDDFT approaches, and proposes a way to 
describe the latter in terms of the former. In Section ITVl 
we derive expressions for steady-state flux, while Sec- 
tion|V]deals with optical response of the current-carrying 
junction. Section |VT1 summarizes our findings. 



II. MODEL 

The model we employ consists of a molecule coupled 
to two contacts (left L and right R). Each contact is 
a reservoir of free charge carriers at its equilibrium, i.e. 
characterized by its own electro-chemical potential /j,k 
(K = L,R). We introduce second quantization field op- 
erators for electrons in the molecule ip a and in the con- 
tacts ip a ,K (c is electron spin and K = L, R) with the 
usual anti-commutation relations 

{^ 1 (r^;R);^(^;R)} = ^ 1)OS «(»'i-^) (1) 
{$<ri,tfi(n);V£ 2 ,ir 2 (r2)} = 5 au(T2 S Ku k 2 S (fi - r 2 ) (2) 

and all other anti-commutators being zero. Note that 
molecular field operators depend parametrically on the 

nuclear configuration R = |i? Q | (Born-Oppcnhcimer 

approximation) . 

Many-body electronic Hamiltonian of the system in the 
second quantization is 

H(R) = H L + H R + H m (R) + Vl(R) + V R (R) (3) 
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where the contact Hamiltonian is given by 

h 2 d 2 



2 m 9r? 



Tpa,K{r k ) (4) 



(here and below K = L,R), coupling between molecule 
and contact K is (no spin-flip transitions) 

V K (R) J dr j ' dr k (5) 

[V^(r; R) V a , K (f, rfejR) $a, K (fk) + H.c. , 
and molecular Hamiltonian is 



2m <9r 



E 



eE(t)r 



^(r;R) (6) 



n - r- 2 \ 



T^a 2 (r2; R) ^(fijR) 



Here a indicates sum over nuclei of the molecule, Z a is a 
charge of the nucleus a, and E(t) is external laser field. 
Note that we intentionally started from the Hamiltonian 
in the second quantization form and used explicit parti- 
tion of the system into contacts and a molecule. Since 
this is a standard NEGF approach, it makes the following 
connection to TDDFT clearer. At the same time, NEGF 
partitioning scheme by Caroli et al^Si and partition- 
less approach by Cini^ implemented in TDDFT-NEGF 
schemes for time-dependent transpor t 74 i 75 ' 76 i 77 ' 78 i 79 in 
the steady-state case, which we are going to consider, 
were shown to be equivalent.— Note that below we will 
consistently use atomic units, i.e. h = e = m = 1, and 
drop R, keeping in mind that all the quantities depend 
on positions of nuclei parametrically. 

Now in the spirit of DFT (and TDDFT) theory we re- 
place the true many-body molecular Hamiltonian ^ by 
fictitious single-particle Kohn-Sham Hamiltonian, which 
in the second quantization form is 



H 



(KS) 
M 



dftp^ (f 



+v* c (r,t) - E(t)r Mr) (7) 



Here 



v ext (r) 



R a 



(8) 



is external potential due to electron-nuclear interaction, 



vf(r,t)= / dn 



\r- ri\ 



(9) 



is electron-electron Coulomb interaction, and 



< c (r,t) 



5A* 



5n a {r,t) 



(10) 



is the exchange-correlation potential. Here A^n] is the 
exchange-correlation action, which should be defined on 
the Keldysh contour^ In Eqs. © and (TIT)]) n a (r, t) is an 
electron density of spin a in position r at time t. 

Below we will represent molecular electron subspace 
in some finite basis set where i is some quantum 

number (or set of quantum numbers) and a is spin. As 
a basis one can use atomic or molecular single-particle 
states (orbitals), or any other convenient basis set. In 
this basis set representation 



i 



(ii) 

f!2) 



where ci a (c| CT ) is electron annihilation (creation) opera- 
tor for state ia. Kinetic energy and potentials in Eqs. ([?])- 
(flU)) become matrices in Hilbert space. Similar basis set 
representation will be used for electrons in contacts. 



III. DENSITY MATRIX FORMULATION OF 
NEGF 

Central quantity of interest within NEGF is an electron 
Green function (GF), defined on the Keldysh contour as 

Gijcr (t, t 1 ) = -i<T c c iCT (r) ct CT (r') > (13) 

where r and r' are contour variables, T c is a contour 
ordering operator, and operators Cj CT are given in the 
Heisenberg representation. Applying standard pertur- 
bation theory and assuming that initial correlations died 
out, one arrives at the Dyson equation on the contour 



OT 



G(t,t / ) = 5(t,t')+ I d7iS(r,Ti)Gf(riy) 

(14) 



or relative to the other variable 



G(r,r') 



d 



H = S(t,t')+ j dri G(r, ti) £(ti,t'). 

(15) 

In (|14| and (fl~5|) . H is Hamiltonian of the part of the 
system under study, while E is self-energy (SE) which 
represents influence of all other parts (and processes). 
Here and below we suppress matrix indices keeping in 
mind that Hamiltonian, GFs, and SEs are matrices in 
the Hilbert space. In the model introduced in Section |TT] 
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H corresponds to Hm, Eq.©, and £ represents contacts 
and coupling to them, Eqs. (g]) and ((5) 

Sjj <T (ri,T 2 ) = 2J ^ifc.g gfco- ( T l > r 2) ^fcj.g ( 16 ) 

K—L,R k£K 

where gka(n,T 2 ) = -i < T c c ka (ri) c^fa) > is GF for 
free electrons in the contacts. Generally £ should in- 
clude also contributions from many-body processes, such 
as electron-electron interaction in Eq.([6|), however in an 
attempt to combine NEGF with TDDFT the last is intro- 
duced through potentials © and (fTU|) in the Kohn-Sham 
Hamiltonian (J7J). 

In what follows we will need lesser projection of (fT4|) 
and {H]) 
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(17) 



+ 00 



dti [T, r (t,t 1 )G < (h,t') + E < (t,t 1 )G a {t 1 ,t')] 



G<(t,t') 



dt' 



(18) 



dh [G< (t, h) Z a (h , t') + G r (t, h) £< fa , <')] 

> 

and retarded projection of (fT4|) 



at 



G r (M') = S(t-t')+ 



dt 1 Tr{t,t 1 )G r {t 1 ,t l ) 
(19) 



while advanced GF is 

G a ija {t,t') = [G r jia (t',t)Y 



(20) 



Note, t and i' in (|17p - (j20|) are time variables. In order to 
get time-dependent solution (without initial correlations) 
within NEGF one has to solve (fTT)) and (fTi?)) simultane- 
ously, while SE should represent also electron-electron 
interaction. 

While in NEGF one deals with retarded and lesser 



Gf ja {t,t')=i<cUt')cio{t) > 



(21) 



GFs, central object of TDDFT is single-electron density 
matrix (DM) 

PiAt) =< 4<r(*) M*) >= - iG ^(t, t) (22) 

Note that rigorous foundations of TDDFT formalism^ 
establish correspondense only between densities of real 
and noninteracting (Kohn-Sham) systems. However 
in practice non-diagonal elements of the Kohn-Sham 
DMs are also used in calculations of optical proper- 
ties. Moreover, approaches explicitly utilizing DM within 



time-dependent functional theory (TDDMFT) were also 
proposed^S^Z. 

Time evolution of DM 



d_ 

dt 



Pit) 



d G <(t,t') + -^G<(t,t') 



dt 



(23) 



can be obtained from Eqs. (11711 and (fT8|) (approximately) 
expressed in terms of G < (t, t) only. In order to do so we 
employ generalized Kadanoff-Baym ansatz (GKBA)^ in 
the right side of (JT7J) and (fT5)) . thus partly loosing non- 
locality in time. This leads to 



4c < (M)-[ff;G < (M)] 
at 



(24) 



/+oo 
dh [Z r (t,t 1 )G<(t 1 ,t 1 )-G<(t 1 ,t 1 )Z a (t 1 ,t)] 
-oo 



-foe 



dh [G r (t, *!)£<(*!,*)-£<(*, h) G a (h,t)] 



Note that in the absence of contacts (i.e. when all SEs 
are zero) and substituting H^f S \ Eq.®, in place of H 
we recover the standard TDDFT formulation in terms 
of DM4i Ea.([24]) together with ([T9j) and ([20j) is the ap- 
proximate formulation for NEGF in terms of DM evolu- 
tion. Below we use these expressions in order to get first 
a steady-state transport through the junction, and then 
optical response of such current carrying junction to an 
external laser field. 

Summarizing, NEGF-TDDFT is superior over NEGF- 
DFT due to ability of treating both time-dependent 
transport and/or optical response of current-carrying 
molecular junctions. Still, as is indicated above, it misses 
nonlocality in time due to GKBA applied and keeps lim- 
itations of TDDFT. Also fundamental question of com- 
bining the two ideologically different schemes remains. 



IV. STEADY-STATE CURRENT 

First we consider steady-state current through the 
junction in the absence of an external field, E(t) = 0. 
In the steady-state situation GFs and SEs in (fl"9]) and 
depend on the time difference only, and DM be- 
comes time- independent, p — —iG < (t = 0). In order 
to simplify notation we further consider contacts within 
wide-band approximation (WBA), when 



^ r (h-t 2 )~--rs(h-t 2 ) 



(25) 



Here 



rv = 27r Y Vik " Vk ^ E - £k ^ = r ^ + r 



R 



K=L,RkeK 



(26) 

is assumed to be constant independent of energy E, e ko 
is energy of the state ka. 
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In this case Eq. (|24|) yields (for brevity we write G < 
keeping in mind that this is G < (t = 0)) 



i^-G<-[H;G<]+ l -{T;G<} = 



dt 

/• + OO 



(27) 



/+oo 
dh [G r (t - ti) E<(ti - t) - E<(t - ti) G a (ti - t)] 
-oo 

where [...;...] is a commutator, while {...;...} is an 
anti-commutator. Note, ij^G < = and is written 
here only in order to keep similarity to the structure to 
Eq.(f2"31). Here 



H, 



/ tjcr 



■xc 

ija 1 w ij(T 



(28) 



with 



dr4>* a {r) 



-v ext {r) 



(ijo-\nmcr')p m na' 



(29) 
(30) 



(ija\nma')= \ dr\ \ df^^in) <j>j a (fi) (31) 

1 



r0n<7'(^2)</>m<7'(^2) 



\ri - r 2 \ 

Note that v xc generally is not a ground state xc func- 
tional. In order to estimate it one can follow approach de- 
scribed in Ref. [92| and utilize formal equivalence (in par- 
ticular cases) between partitionless TDDFT and parti- 
tioned NEGF schemes for the case of established steady- 
state current^ For details of approximate way to esti- 
mate v xc see Appendix 1X1 

Fourier transform (FT) of lesser SE entering Eq. (|2"7| 



is 



X<(E) = i[f L (E)T L +f R (E)T R ] 



(32) 



where fx(E) is the Fermi distribution in the contact K, 
and retarded GF G r can be obtained from FT of (fT9l) 



G r {E) = [E-H + iT/2] 1 



(33) 



Note that chemical potentials in the contacts should be 
shifted to take into account xc response to bias induced 
fields 

Integral version of Eq. |27|) (Keldysh equation) 

/+oo irp 
— G r {E)Y, < {E)G a {E) 1 (34) 
-oo 27T 

together with Eq.® yields the standard NEGF-DFT 
approach to steady-state transport. Note however that 
in the last case v xc is substituted by ground state xc 
potential and xc corrections to chemical potentials of the 
contacts are neglected. 



V. LINEAR OPTICAL RESPONSE 

Here we consider linear optical response to weak ex- 
ternal laser field E(t). In contrast to previous TDDFT 
considerations^^ the response is calculated on top of 
non-equilibrium steady-state (rather than ground state) 
of the system. 

Upon introduction of the time-dependent external field 
Eqs. (HU and © yield (within WBA) 

i -G < (<, t) - [H; G<(t, t)\ + - {r ; G<(t, t)} = (35) 

-foe 

dh [G r (Mi)£<(ii -t) -Z < (t~t 1 )G a (t 1 ,t)] 

-oo 

G r {t,t')=G r {t-t') (36) 

/+oo 
db x G r {t-t 1 ) [H-H] G r (t u t') 
-oo 

Note, in Eq. (J3SI) wc assume that contacts are not in- 
fluenced by an external field (SEs are the same), and 
Eq. (f3"6"|) is an integral variant of Eq. (TT9"|) . Here 



Hi 



H 



(KS) 
M 



hija + vfUt) + vfUt) - fl ija E{t) 



(37) 



with hiju defined in Eq. (|29|) . vfj a (t) defined similar to 
Eq. ([30|) with p replaced by p(t), vff a (t) is xc potential 
when both bias and external optical field are applied to 
the junction, and the matrix element of molecular dipole 



(Uj<T= / dr(j>* ia (f)r(j> ja {r) 



(38) 



H is defined in Eq. 

Assuming external optical field is a weak perturbation, 
we linearize Eqs. (|3"5"j) and (|3"6"|) in E(t). Introducing re- 
sponse to the field 

6G<(t,t) =G < (t,t)-G < = iSp(t) (39) 
5G r ' a (t, t') = G r < a (t, t') - G r - a (t - t'), (40) 

linearizing Eqs. ([55)) and (f3"6")) in the field and response to 
it, and subtracting Eq. (|27p from the linearized version of 
Eq. ([35l) we obtain 



ij t 5G< (t, t) - [H; 8G< (t, t)\ + t -{T; SG< (t, t) } 



[Sv cl (t) + 8v xc (t); G<] + E{t)[p- G<] = (41) 



+00 



dh [5G r {t,t 1 )Z < {t 1 -t)-Y< < (t-t 1 )5G a {t 1 ,t)] 

■oo 

/+oo 
db x G r {t-t 1 ) (42) 
-OO 

x \6v cl (h) + Sv xc (h) - E(ti) p] G r (h - t') 
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Here 



8vf j<r {t) = v? i<r (jk)-v. 



(ij(x\nma')8p mn(T >(t) 
(43) 



E 



+ 00 



^1 fijcr,nmcr' (* — £ l) <>Pmncr> (h) 

(44) 



where superscript "r" on xc response function / r (£ — ii) 
indicates its retarded character, i.e. for t < ti f r (t—ti) = 
0. Since steady state can be (approximately) mapped to 
an effective equilibrium , 89 ' 90 ' 91 Sv xc (t) in principle can 
be estimated following procedure of Ref. [92] when equi- 
librium is understood as an effective equilibrium corre- 
sponding to a particular steady-state. 

Substituting Eq. (|42|) and its advanced analog into 
Eq. (|41|) and rearranging terms so that all source terms 
(terms containing external field) arc in the right-hand 
side, we obtain final result, which in Liouville space is 
given by 



9 T ' 

i— - L Q 



at 



S P <(t)- 



-foe 



dh L{(t-ti)5p(ti) 
dtt&it-t^Efa) (45) 



where superscript "r" again indicates retarded charac- 
ter of the corresponding functions. Fourier transform of 
Eq.(|i5]) is 



(46) 



[u-Lq- L{(lu)} Sp(cu) = £ r (uj) E(w) 
Here Markovian part of the Liouvillian is 



[L . 



% j cr, mncr 



(47) 



H, 



njcr 



2 n 3 a 



frequency-dependent Liouvillian is 

^ ] \_f ipir ,nmcr' (^) Ppjer Piper fpjcr.nmcr' M] 



*E 



dE 



{G r ipa (E + uo)G r qra (E)^< a (E) (48) 

^irA E ) G rqa( E ) G pja( E ~ w ) 

x [(qpa\nma') + f qp ^ nm( ,, (w)] } , 



and frequency-dependent coefficient of the source term is 

^ij(f^) = ^ ^ \Pipo Ppja ~ Piper Ppjcr] 



v , r°° dE 



[G\ p(T {E + u,) G r qra (E) E<. ff (S) p pqa 

^frA E ) G rqa( E ) G pia( E ~ W ) Pqpo\ 



(49) 



Since in the linear optical response only particle-hole, 
Pno > Pjjcr, and hole-particle, pu a < pjj a , elements of 
Spij a (t) are nonzero,— size of the matrix equation (|46p 
can be reduced. Following Ref.[93| we order the basis such 
that i < j for pu a > pjja, and divide Sp(t) into particle- 
hole and hole-particle parts. Introducing matrices in the 
ordered basis, so that i < j and m < n 



Aija,rnncr'{uj) - [L + £l(w)] ijV>r)m(T , 
Bija,rnncr'{uj) = [L + L{ (w)] ijV>nm(T / 

and using general relations 

5pij a (u) = Sp* ia (-u) 
r r l _ _ it 1* 

L ^Hja,mna' L ®\jio,nm<7 ! 

[ L l{ W )]ija,mna' = ~ 1^1 *jia,nma> 



we get matrix equation of the form 

cj - A(w) — B(w) <fy>(<^) 
B*(-w) w + A*(-w) J [5p*(-uj) 



(50) 
(51) 

(52) 
(53) 
(54) 

(55) 



(56) 



where 



S r (w)=£ r (u})E(u) 



(57) 



Eg. ([56]) is reduced dimension form of Eg . (|4"6"|) . In the case 
when contacts are absent and real molecular orbitals are 
chosen as a basis, matrices M = {A, B, S r } satisfy 



M*(-lu) = M{uj) 



(58) 



and Eq. ([56"]) reduces to the result derived within TDDFT 
approach previously^ 

Polarization of the molecule is 



P a (t) =Tr[p a Sp(t)} 



E 





dtxFQit-t^E^tx) 
(59) 

where a, /3 = (x,y,z), R^lit — t\) is the time-domain 
linear response function (tensor in x,y,z), and trace is 
over molecular electronic basis (TV [. . .] = YliA- • •]<»)■ 
Eg. ([46[) yields for the frequency-domain linear polariz- 



ability (FT of (t - h)) 



* a p(uj) = -Tr |/.t Q [u - L Q - L\(lo)] 1 £^{u)^ (60) 
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Excitation energies (continuous spectrum of excitation 
energies) and oscillator strengths can be obtained from 
the poles and residues of the polarizability. In the ab- 
sence of contacts Eg. (fBDj) reduces to the result derived 
previously^ 



VI. CONCLUSION 

We discuss correspondence between NEGF and 
TDDFT approaches, and derive approximate DM rep- 
resentation for NEGF equations. Thus derived NEGF- 
TDDFT approach is superior over NEGF-DFT due to 
ability to treat both time-dependent transport and/or 
optical response of current-carrying molecular junctions. 
However it partially misses nonlocality in time, inherent 
to NEGF. Also fundamental question of combining the 
two ideologically different schemes (NEGF and TDDFT) 
remains open. We further propose a practical scheme for 
calculation of linear optical response of current carrying 
molecular junctions. Situation considered corresponds to 
the case when an electron transport through the junction 
is much faster than the characteristic time of an external 
laser field (e.g. a pulse length). 

First we derive expression for steady-state cur- 
rent through the junction within NEGF-TDDFT 
scheme and discuss its correspondence to the TDDFT- 
NEGF approaches to time-dependent transport of 
RefsM£&&ZLZ&22.. Formal equivalence of NEGF- 
TDDFT and TDDFT-NEGF schemes in the case of es- 
tablished steady-state flux permits to propose a way to 
estimate xc potential, and thus go beyond ground-state 
xc potential implemented in the standard NEGF-DFT 
transport schemes. 

After that we consider response of such current car- 
rying junction to weak external laser fields. We derive 
expressions for linear response and polarizability, which 
in the absence of contacts reduce to previously obtained 
TDDFT results for optical response of isolated molecules. 
Presense of the contacts introduces memory in both Liou- 
villian and source term due to external field. Equation for 
response DM is then expressed in reduced (particle-hole 
and hole-particle) basis, which allows direct comparison 
to corresponding isolated molecule expression obtained 
previously. 

The proposed approach allowing to calculate optical 
response of current carrying junctions is a first step in 
direction of ab initio calculations of such kind, with the 
final goal to go beyond model based studies of optical 
properties of current-carrying junctions available in the 
literature. Development of such ab initio schemes is es- 
pecially important in light of experimental data on opto- 
electronic properties of molecular junctions which start 
to appear. In contrast to optical response of an isolated 
molecule presence of contacts will smear discrete excita- 
tion spectrum into continuous one, which is in complete 
analogy to smearing of discrete energy spectrum of iso- 
lated molecule into continuous density of states upon at- 



taching to the contacts. Presence of contacts also allows 
for charge transfer transitions upon optical excitation - 
process absent in the isolated molecule case, and thus 
changing the optical resonse of the molecule. Finally, 
nonequilibrium character of the junction will alternate 
optical response. An obvious change is presence of elec- 
tronic density in an extended energy region, defined by 
difference of electro-chemical potentials of the two con- 
tacts, opposed to well-defined ground state energy in the 
case of isolated molecule. For example, this may lead 
to appearance of additional (inverse) Raman scattering 
channel (situation when initial and final electronic state 
is LUMO rather than HOMO in the isolated molecule 
case), and to interference between the two (normal and 
inverse) channels, as is discussed in detail elsewhere^ 

Consideration of situation, when time of electron tun- 
neling and characteristic time of external field change are 
comparable, is a subject of future studies. 
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APPENDIX A: APPROXIMATE WAY TO 
ESTIMATE v xc 

Here we discuss an approximate way to estimate v xc for 
current-carrying molecular junction in steady-state. We 
follow approach by van Leeuwen^, where connection be- 
tween xc potential and standard diagrammatic technique 
was established. Ref. [92| considers a system both within 
full many-body approach and within TDDFT, where true 
Hamiltonian H is replaced by fictitious single-particle 
Kohn-Sham Hamiltonian H^ KS ^ with some xc potential 
v xc to be found. Then adiabatic switching on of the in- 
teraction, H — H^ KS \ is considered, which leads to con- 
nection between full many-body GF G and single-particle 
Kohn-Sham GF Gks m the form 

G(fi,Ti;r 2 ,T 2 ) = G K s(ri,T 1 ;f2,T 2 ) 
+ j " dr 3 J dr 3 J dr A J dT4,G K s(ri,Ti;r 3 ,Ta) (Al) 

^xc(r3,T 3 ;f 4 ,T 4 ) - S(f 3 - r 4 )<5(T 3 ,T 4 )i; (:,;c) (f3,T3) 

x G(r i ,T i ;r 2 ,T 2 ) 

with T. xc being xc part of the SE due to electron-electron 
interaction. Since the two approaches are assumed to give 
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the same electron density 

n(r, t) = -iG < {r 1 t; r, t) = -iG< s {r 7 t; r, t) (A2) 



lesser projection of (jAlj) provides a (self-consistent) ex- 
pression for xc potential v^ xc ^ in terms of SE E xc 



dh [G KS {r,t;r u h) v xc (n , h) G(f u h; f, t)} < 

(A3) 



— oo 

/+oo /> /"+00 

di\ \ dr2 / dt2 
-oo J J— oo 

[Gks (r , t ; r i , ti ) (r i , t i ; f 2 , t 2 ) G(r 2 , t 2 ; f, i )] < 
where 

[G^s w ( " c) G]< = G ks v {xc) G a + G r KS v (xc) G< (A4) 
[G KS Ere G] < = Gf iS E" c G a + G r KS E< c G a ( A5) 
+ 5 S^ c G < 



Eg. (|A3[) is our starting point. 

Now we consider situation of established steady-state 
current in the system. In this case it was shown that 
v xc becomes constant far from the device region.— Ac- 
cordingly, we formally partition the system into 3 parts: 
two contacts (A" = L,R), where v xc is constant, and 
device (molecule, M), where it varies in space. Then in- 
tegral over space in the left of Eq. (|A3[) splits into two 
parts (integration over contacts Vk and over device Vm), 
while double space integral in the right will have 4 con- 
tributions: two over Vk and Vm, respectively, and two 
mixed ones. The last may be neglected (especially when 
device is well separated from the contacts) taking into 
account relatively short range of xc interaction. This 
approximation yields additive (in K and M) structure 



of Eq. (|A3[) . and obviously one can equate corresponding 
parts separately. As a result one gets for the device re- 
gion only expression similar to Eq. (|A3[) . where v^ xc ^ is 
time-independent and where GFs are restricted to the 
device region only. As was demonstrated in Ref. [zl the 
last are given formally from the usual NEGF expressions 
for GFs of the device region in terms of SEs due to cou- 
pling to the contacts, when appropriate constant shift of 
chemical potentials in the contacts due to xc response is 
taken into account. Moreover, while generally this shift 
may depend on the history of appearance of the steady- 
state, in LDA it depends only on the instantaneous local 
density and has no memory at all^ i.e. instantaneous 
(and constant in time) potential and instantaneous (and 
constant in time) density uniquely define each other. In 
this case we can assume that the correct (shifted by xc 
response) chemical potentials in the contacts are some 
predefined boundary conditions, as is usually done within 
NEGF, which determine situation in the device region. In 
other words we assume equivalence of partitioned NEGF 
and partitionless TDDFT approaches in this case. Then, 
after introducing some basis within the device region (so 
that GFs, SEs, and v xc are matrices in Hilbert space), 
Eq-CMI yields 



/ [G KS (E)v xc G(E)} < = 
^ [G KS (E)^ XC (E)G(E)} < 



(A6) 



with GFs given by usual NEGF expressions. Eq. (|A6|) 
is an approximate way to estimate v xc in the steady- 
state situation. Explicit expression for v xc is obtained 
by rewriting Eq. (|A6|) in Liouville space. 
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